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Abstract 

By using the three-term recurrence equation satisfied by a family of orthogonal poly- 
nomials, the Christoffel-Darboux-type bilinear generating function and their asymptotic 
expressions, we obtain quadrature formulas for integral transforms generated by the clas- 
sical orthogonal polynomials. These integral transforms, related to the so-called Poisson 
integrals, correspond to a modified Fourier Transform in the case of the Hermite poly- 
nomials, a Bessel Transform in the case of the Laguerre polynomials and to an Appell 
Transform in the case of the Jacobi polynomials. 
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1 Introduction 



Certain integral transforms with a Christoffel-Darboux-type bilinear generating function 
for a classical orthogonal polynomial in the kernel, termed Poisson integrals [H[2], appear 
in the expansion of functions in terms of the classical orthogonal polynomials and, un- 
der a suitable change of variable, they become well-known integral transforms for certain 
limit values of the expansion parameter. In this sense, Mehler's formula yields a modified 
Fourier transform in the case of the Hermite polynomials, the Hille-Hardy formula gives 
a Bessel-Hankel transform in the case of the Laguerre polynomials and Bailey's bilinear 
generating function produces an Appell transform in the case of the Jacobi polynomials. 
We show in this paper that each of these integral transforms has a quadrature formula of 
gaussian type. This was done by following mutatis mutandis the approach of [U [5], where 
a quadrature formula for the Fourier transform and other one for the Hankel transform 
were obtained by using the recurrence equation for the (generalized) Hermite polynomi- 
als and their zeros as well as the Christoffel-Darboux formula, Mehler's formula and the 
Hille-Hardy formula. The generalization of this procedure provides a method to yield at 
one time new integral transforms and their quadrature formulas. In the case of a family of 
orthogonal polynomials & = {Pk(x), k — 0, 1, . . .} this method requires a suitable bilinear 
generating function of & as well as some knowledge about the asymptotic behavior of 
the polynomials and their zeros. The kernel of the integral transform becomes related to 
the generating function. We apply this technique to the classical orthogonal polynomials 
to obtain quadrature formulas for their corresponding Poisson integrals. Thus, in Sec. 
[3] we present a quadrature for a modified Fourier transform generated by the Hermite 
polynomials, a quadrature for a modified Bessel transform yielded by the Laguerre poly- 
nomials in Sec. HI and finally, in Sec. [5]we introduce a quadrature for an Appell transform 
generated by the Jacobi polynomials. This integral transform has also been obtained in 
0- 



2 Outline of the method 

Let = {Pfc(x), k — 0, 1, . . .} be one of the families of classical orthogonal polynomials 
(Hermite, Laguerre or Jacobi) satisfying the recurrence equation 

A n P n+1 (x) + B n P n (x) + C n - 1 P n - 1 (x) = xP n (x), n = 0,1,2,..., (1) 

where P-i(x) = 0. Let k n denote the coefficient of x n in P n (x). As it is well-known [9l[T0]. 
from ([I]) follows the Christoffel-Darboux formula 

V slP n (x)P n (y) = kN -f N -' P ^)Ps-i(y) - PN-i(x)P N (y) ^ (2) 

n= kN x ~ y 
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where s n is determined by the reciprocal of the norm 

\ 1/2 

P 2 (x)duj(x)) (3) 

up to a numerical constant independent of n. Here, (a, b) is the orthogonality interval of 
and du(x) is the corresponding non- negative measure. The recurrence equation ([T]) 
can be written as the eigenvalue problem M^P = xP, where (Moo) nfc = A n _i5 n+ i <k + 
B n -i8 n> k + C n ^n,k+ii n, k = 1,2, . . ., and P is the vector whose nth entry is P n _i(x). Let 
us now consider the eigenproblem associated to the principal submatrix of dimension iV 
of Moo. This part is a well-known technique [6j [7J E] to yield gaussian quadratures: the 
three-term recurrence equation is rewritten in matrix form to obtain orthonormal vectors 
of R whose entries are given in terms of the values of P k (x), k = 0, 1, . . . , N — 1, at the 
zeros of Pn{x). To proceed, we take a similarity transformation to symmetrize M^. The 
diagonal matrix S = diagjso, si, • . . , } whose elements s n are given by generates a 
symmetric matrix SM^S^ 1 whose principal submatrix of order N, denoted by M, has 
elements given by 



M nk = \j An-xCn-xbn+xfc + -Bn-l^n./e + \j A n ^ 2 C n -2Sn t k+l, U,k = 1,2, N . 

The recurrence equation ([1]) and formula ([2]) can be used to solve the eigenproblem 

Mu k = x k u k , k = 1,2, . . . , N. 

The eigenvalues x k are the zeros of Pn{x) and the nth entry of the kth eigenvector u k is 
given by c k s n -iP n -i(x k ), n = 1,2, . . . , N, where c k is a normalization constant which is 
obtained from (j2j). Thus, the orthonormal vectors u k , k = 1,2, . . . , N, have components 

viW = H p Pn 'f x t ( m ' » = i.—*- (4) 

\ s iV-l/ ^/kN-iPN-^XkjP^Xkj/kN 

Note that the product k^^iP^^i{x k )P' N {x k )/kjq is always positive since it is 1/ s N _]_ times 
the squared norm of the vector 

(s P (x k ), SiPi(x k ), • • • , SAT_iPiv_i(x fc )) . 

Let U be the orthogonal matrix whose kth column is u k and D be the diagonal matrix 
D = diagjl, z, z 2 , . . . , z N ~ x }, where z is a complex number. Then, the matrix T = U~ l DU 
whose elements are explicitly given by 

N-l 



T i k = znu n{Xj)u n {x k ) 



n=0 

\k^ /k N _i\ 



s 2 N -i\/\Pn-i (x^P^x^Pn-! {x k )P' N {x k )\ f 



^4P n (x,)P„(x fe )^, (5) 
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is the Discrete Transform associated to the corresponding Integral Transform. To show 
this, we use some asymptotic properties of the zeros of the classical orthogonal polynomials 
shown for each case in the next sections. First, we note that the asymptotic expressions 
for P/v_i(x) and P' N (x) in the oscillatory region, evaluated at the zeros of Pn{x) satisfy 

the formula 

P N ^{x k )P' N {x k )\ = + 0(N-»). (6) 

The function g(x) is positive in (a, b) and it is related to the weight function of . The 
constant > is independent of k and the exponent /i can be taken positive. The zeros 
of Pn(x) contained in the fixed interval (c,d), a < c < d < b, become evenly spaced for 
large values of N under a bijective mapping y = a(x), i.e., 

Ay k = Aa(x k ) = a(x k+l ) - a(x k ) = X N + O(iV^), (7) 

where \n does not depend on k. Furthermore, we have that 

\k N /k N _i\ 



<? 2 a 2 
b N-l a N 



AAa(x k )+0(N-»), (8) 



where A is a numerical constant independent of N and k. Since the Christoffel-Darboux- 
type bilinear generating function G(x, y, z) of the set & satisfies 

oo 

J24Pn(x)P n (y)z n = G(x,y,z), (9) 

n=0 

for (x, y) G (a, b) x (a, b) and z in a suitable domain T> of the complex plane, the element 
Tj k [cf. Eq. (jSJ)] has the limiting form 

Ag(x j )g{x k )G{x j , x k , z)Aa(x k ), (10) 

for large values of N. Therefore, 

N N 

^T jfc /(x fc ) -> A^2g(x j )g(xk)G(x j ,x k ,z)f(x k )Aa(x k ), N -> oo. (11) 

k=l k=l 

The sum of the right-hand side of this equation is the Riemann-Stieltjes sum of f(x) with 
respect to <j{x) and tends to the integral transform 

T[/(x); y, z] = f K(x, y, z)f{x)da(x), (12) 

J a 

where K(x,y,z) = Ag(x)g(y)G(x,y, z). Therefore, Eq. (iTTj) becomes the quadrature 
formula 

,6 



/ K(x,yj,z)f(x)da(x) ~ ^T jk f(x k ), (13) 
Ja k=i 



for 2; G T>. 
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3 A modified Fourier transform 



As it is known, in the case in which P n (x) is the nth Hermite polynomial H n (x), the 
parameters k n and s 2 are 2 n and l/(2 n n!) respectively and the nodes Xk are zeros of 
H N (x). Taking into account that \H' N (x k )\ = (— l) N+k Hjtf-i(xk), the nth component of 
the kth orthonormal vector u k given by Eq. (j3J) becomes 

^^.^(^^Y 2 ^, B=1 ,..., w . (14) 



jV(n-l)! y #Ar-i(z fe 
Therefore, the elements of the Discrete Transform T, denoted in this case by Tjj, are 

9 N ~ X (N — 1 V N ~ 1 7 n 

(T H ) jh = (-iy +k NR ( [ m \ S £ ^#n(^)#n(^)- (15) 

NH N -. 1 (xj)H N -i{x k ) ^ 2"n! 
By using the asymptotic formula [9] 

= ^^y ^ /2 (cos(V2iVTl x - A%/2) + ^(iV- 1 / 2 )) , 

we obtain that g(x) = e~ x2/2 , a(x) = x, Ax k = n/V2N, A = 1 [cf. Eqs ©-©] and the 
expression for G(x, y, z) is determined by Mehler's formula. Thus, we get the quadrature 
formula 

N 



/oo _ _ 

K H (x, yj, z)f(x)dx ~ ^2(T H ) jk f(x k ), (16) 
■°° fe=i 



where 



1 / (l + z 2 )(x 2 + 2/ 2 )-4x?/^ 
K H [x,y,z) = — = exp ; ^ , z < 1. 

1 ; ^(i-^) p v 2(1 -z 2 ) ; 11 

The argument of the exponential is imaginary if z ^ ±1 is on the unit circle. In partic- 
ular, if z = ±z, Kn(x,y,z) = exp(±ixy) / y/2n , and (fTET) becomes a quadrature for the 
Fourier Transform. This formula has been previously obtained in [4] were some numerical 
examples also has been given. 



3.1 An example 

We fix the point yj where ( 1T6I) is evaluated and thus, the integral becomes a function only 
of z. For z real, the kernel is a real exponential but even so the Fourier transform can be 
recovered by taking f(x) as a complex exponential. In this example we take yj = (N 
odd) and f(x) = exp(-ix). Thus, ([TBI) becomes 

O _ 2 N 

e ~W^ ~ ^(T^^e-^. (17) 
fc=i 
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In Fig. 13.11 we show the numerical computation of this formula. 




4 A Bessel transform 

We take now the Laguerre polynomial Ln(x) as P n (x). Thus, k n = (— l) n /n\ } 
n\/T(a+n+l) and the nodes x k are the N zeros of L^\x). We have now that \L £ (xk)\ 
(— (N + a) I 'x]JJ^_ x {x^) and the vector u k given by (TJ]) becomes 

u Ax k ) = (-1)*+* f (»- l)»r(a + JV)*» y* L&M 



> l) j+k T(a + N)JxIxZ riz n 



The matrix T, denoted now by Tl has components 

N-l 

iV!(iV + aji^)^) ~ r(a + n + 1 
The use of the asymptotic formula [9] 

jya/2-l/4 e x/2 



0p x a/2+l/4 



cos(2v / iV^ - «tt/2 - tt/4) + (iVx)~ 1/2 C>(l)) 



yields #(x) = a^+^e - */ 3 , a ^ = ^ Aa(x k ) = n/2y/N, A = 1/2. The generating 
function G(x, y, z) is given by the Hille-Hardy formula. Thus, we get the quadrature 
formula 

poo x 

/ K L {x, yj ,z)f{x)dx ~ V(T L ) ifc /(x fc ), (20) 
■ / ° fc=i 



where 

-a/2 



Izl < 1, 



and is the modified Bessel function of the first kind. Note that we have used the 

fact that da(x) = x~ 1 ^ 2 dx/2. 

For z — — 1, fl20]) becomes a quadrature for the Hankel transform 

2 ^ I x ~ 1/4j <*(Vw)f(x)dx - ^(T L ) jfc /(x fc ), (21) 



k=l 



which is written in an nonstandard way. 



4.1 An example 



The integral 



" ' xJ a {cx)I a {bx)dx 



e 4a J ( 

2a v 2a' 



(22) 



is given usually as a Hankel transform [13J, but it will be used here to test (1201) . To this 
end, we make the change of variable x — > x 2 in the left-hand side of (1201) and take 

1 + 2 



2(1 



/(x 2 ) = ^fxJ a (cx). 



Thus, ( 120]) becomes 

2z-/ 2 1/4 

-j— y/ exp 

l + z J 



1 + 2 



c + Vi 



l + z 



N 



J,, ( \ ^ y^(T L ) jfc Xfc /4 J a (c v /a^). 



fc=i 



The numerical output of this quadrature is illustrated in Fig. 14.11 Note that 
exact formula if z — 1 since T becomes the identity matrix for this value of z. 




o 

X 0.05 
W 



15 20 

y 



25 30 35 



(23) 



is an 



Figure 2: The left-hand side of ( l23l) is plotted versus the right-hand side. The values of 
the parameters are N = 30, a = 0, c = 2, and 2 = 1/10. The norm of the error is 0.0036. 
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5 An Appel transform 



Proceeding as before, we take now the Jacobi polynomial P^'^\x) as P n (x). Therefore 



{2n + a + (3)\ 
2 n (n + a + (3)nV 



(2n + a + (3 + l)n\T(n + a + (3 + 1) 



T(n + a + l)r(n + (3 + 1) 



In this case, the nodes x k are the zeros of P]^ (x). We have now that 



y(«,/9) 

N 



M\ 



1)N+k 2(N + a)(N + (3) p W0 M 



(2N + a + (3)(l-xl)' 
and the components of the kth orthonormal vector u k are 

(n - l)!(2n + a + (3 - 1)(2N + a + (3) 2 T(N + a)T(N + (3) 



U n -l(Xk) 



-1) 



N+k 



AN\(N + a)(N + (3)T(N + a + (3 + l)T(n + a) 



T(n + a + (3) 



1/2 



X 



n 



1, . . . , N. 



r(n + /?) 

Therefore, the components of the Discrete Transform T, denoted here by Tj are 



(24) 



(Tj) 



JV-l 



fc (2N + a + (3) 2 T(N + a)T(N + (3) yji 1 ~ tfK 1 ~ x l) 
AN\(N + a)(N + (3)T(N + a + (3 + 1) P^( Xj )P^(x k ) 



x 



w!(2n + a + (3 + 1)I> + a + (3 + 1) ^p^,/?)^.^,/?)^^ 



n=0 



T(n + a + l)T(n + (3 + 1) 
The use of the asymptotic formula [9] 

i cos [(N + z&±)e -(« + !)!] 



(25) 



Pi a,/3) (cos, 



7T 



(sm |) (cos 2) 



0(AT- 3 / 2 ) 



yields g(x) = (1 -x) a / 2+1 / 4 (l + x) /3/2+1/4 , a(x) = arccos(:r), Aa(x k ) = n/N, A = 2~ a ~ 13 - 1 . 
The generating function G(x, y, z) is given by Bailey's formula [TT] 



Gj(x,y,z) 



x 



r(q + /? + 2) 
r(a + 1)T(P+ 1) (l + z) a +P+ 2 
a + (3 + 2 a + (3 + 3 



;a+l,/3 + l; 



z(l-x)(l-y) z(l + x)(l+y) 



(1+*) 



where F 4 (a, 6; c, d; £, 77) is the fourth Appel's hypergeometric function of two variables 
If r and s are the radii of convergence in £ and 77 respectively, then r 1 / 2 + s 1 / 2 = 1 and x, 
y and z should satisfy 



z{l - x)(l - y) 
(1 + z) 2 



< r, 



z(l + x)(l + y) 
(1 + z) 2 



< s, 
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in order to have the quadrature formula 

•1 N 



[ Kj(x, yj ,z)f(x)dx~J2(Tj) jk f(x k ), (26) 
J ~ l k=i 



where 



j\j[x,y,z)- 2 a +' 3 +W l-x 2 <^j{x,y,z). 

Note that we have written the integral in terms of dx instead of da(x). 



5.1 An example 

It is shown in [3] that0 

J (1 - x) a (l + xfP^\x)F i {a, ft a + 1, (3 + 1; £, rj)dx 

_ 2^I> + 1W +!)(! + g , p( ^M f27) 
(a + /3 + l)r(a + ^ + l) (1-*) " 1 J 

where 

_ _ a + P + 2 - _ a + g + 3 z(l - g)(l - y) _ z(l + g)(l + y) 

"~ 2 ' ^~ 2 ' 4 ~ (1 + ^) 2 ' // ~ (1 + ^) 2 ' 

By choosing a suitable integrand, our quadrature formula (T26]) becomes 

AT 

^(l-yjO^^Cl+yjO^^P^^i) ^ 5Z(T J ) ifc (l-x fc ) a/2+1/4 (l+x fc ) /3/2+1/4 Pi a ' /3) (x fc ) 

fc=i 

(28) 

for this case. The numerical output of this quadrature is illustrated in Fig. 15.11 This 
equation is an exact formula if z — 1. 



1 We have corrected a misprint in this formula. 
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xlO- 4 




y 

Figure 3: Plot of the left-hand side against the right-hand side of (128]) for N = 50, 
a = (3 = 0, z = 1/4 and n = 5. The norm of the error is 0.00011 in this case. 

6 Final Remarks 

There are important facts supporting this method for finding quadrature formulas for the 
integral transforms (Poisson integrals) associated to a family of orthogonal polynomials. 
These facts are given by formulas (EI)-([9]) . An important one, is that the kernel of the inte- 
gral transform is determined by the product of the bilinear generating function G(x, y, z) 
and g(x)g(y). The bilinearity of the kernel explains why the associated quadrature for- 
mula can yield exact results if the function to be transformed is chosen properly in terms 
of elements of the orthogonal family. 
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